Robust Standard Errors

1 High Level Summary

We obtain different estimates for the standard deviations of the parameter estimates. In particular, we compare estimates under the assumption of correct model specification with robust ones.

Under correct specification, the information matrix matrix may be calculated as the

• outer product gradient (OPG) of the scores
• the Hessian obtained as output for stats::optim()
• the analytic version of the Hessian obtained from numDeriv::hessian()

Allowing for misspecification, the robust standard errors are obtained from the covariance matrix \(\Omega^{-1}\left(\hat{\theta_T}\right) I\left(\hat{\theta_T}\right) \Omega^{-1}\left(\hat{\theta_T}\right)\). Again, the Hessian may be calculated as output of stats::optim() or in an analytic way from numDeriv::hessian().

First, we investigate the OPG and the different versions of the Hessians quantitatively and by plotting the values, see section OPG and Hessian and in particular Plot Parameters. Subsequently, we compare the diagonal elements of the OPG and the different versions of the Hessian in [Diagonal Elements]. Next, we evaluate the implied standard deviations for the system and noise parameters quantitatively and with plots in Compare Standard Deviations. Finally, we provide the results, i.e. each part of the parameter estimates together with their robust and non-robust standard errors, see Results.

Some distribution parameters are not estimated very precisely as one can see from Noise Parameters. The sub-matrix pertaining to the distribution parameters is close to singular.

2 Load data

2.1 Data Set

Load Blanchard and Quah (1989) dataset and its dimensions.

DATASET = readRDS(params$FILEPATH_DATASET) %>% 
  as.matrix()
DIM_OUT = dim(DATASET)[2]
N_OBS = dim(DATASET)[1]

2.2 Optimal model

Load the template and the deep parameters for the best model.

PARAMS = read_rds(glue("{params$DIRPATH_MOD}bf_params_12_2.rds"))
TMPL = read_rds(glue("{params$DIRPATH_MOD}bf_tmpl_12_2.rds"))

Fill the template

params_opt_mod = fill_tmpl_whf_rev(theta = PARAMS,
                                   tmpl = TMPL)

Print AR parameters in wide matrix format.

params_opt_mod$polm_ar %>% matrix(nrow = 2)
##      [,1] [,2]       [,3]       [,4]
## [1,]    1    0 -0.1073603 -0.2770627
## [2,]    0    1  0.4921382 -1.0125875

Print MA backwards parameters in wide matrix format.

params_opt_mod$polm_ma_bwd %>% matrix(nrow = 2)
##      [,1] [,2]      [,3]       [,4]
## [1,]    1    0 0.0650218 0.74866119
## [2,]    0    1 0.2593842 0.00153858

Print MA forwards parameters in wide matrix format.

params_opt_mod$polm_ma_fwd %>% matrix(nrow = 2)
##      [,1] [,2]       [,3]       [,4]
## [1,]    1    0  0.5610642  1.6709962
## [2,]    0    1 -0.2414728 -0.2956298

Print static matrix B

params_opt_mod$B %>% matrix(nrow = 2)
##            [,1]      [,2]
## [1,]  0.9357336 0.2888584
## [2,] -0.2234818 0.1819431

Print distribution parameters, \(\lambda, p, q\) for each variable (in columns)

params_opt_mod$distr %>% matrix(ncol = 2)
##               [,1]      [,2]
## [1,]    -0.5243178 0.1492405
## [2,]     1.4924937 1.9211332
## [3,] 85616.3114036 7.8907217

Number of deep parameters parameters

TMPL$ar$n_par
## [1] 4
TMPL$ma_bwd$n_par
## [1] 4
TMPL$ma_fwd$n_par
## [1] 4
TMPL$B$n_par
## [1] 4
TMPL$distr$n_par
## [1] 6

2.3 Indices of System and Noise Parameters

ix_params_sys = 1:12
ix_params_sys_ar = 1:4
ix_params_sys_ma_bwd = 5:8
ix_params_sys_ma_fwd = 9:12

ix_params_noise = 13:22
ix_params_noise_b = 13:16
ix_params_noise_sgt = 17:22

2.4 Likelihood function

Obtain likelihood function for given dataset and template.

ll_fun_sgt = ll_whf_factory(data_wide = t(DATASET), 
                            tmpl = TMPL, shock_distr = "sgt", 
                            use_cpp = TRUE)

Recalculate the optimisation with argument hessian = TRUE.

N_PARAMS = length(PARAMS)

out = optim(par = PARAMS, 
            fn = ll_fun_sgt, 
            method = "L-BFGS-B", 
            lower = c(rep(-Inf, N_PARAMS), rep(c(-1, 1.5, 2), DIM_OUT)),
            upper = c(rep(Inf, N_PARAMS), rep(c(1, Inf, Inf), DIM_OUT)),
            control = list(maxit = 100), hessian = TRUE)

Check output

out$par - PARAMS
##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

3 OPG and Hessian

3.1 Obtain Different Versions

Calculate OPG

ll_grad = ll_whf_factory_grad_helper(data_wide = t(DATASET),
                                     tmpl = TMPL)

J = jacobian(func = ll_grad,
             x = PARAMS)
opg = J/sqrt(N_OBS)
opg = crossprod(opg)

if(params$ZERO_RESTR){
  opg[ix_params_sys, ix_params_noise] = 0
  opg[ix_params_noise, ix_params_sys] = 0
}

opg_inv = ginv(opg)

Calculate Analytic Hessian with numDeriv package.

ha = hessian(func = ll_fun_sgt,
             x = PARAMS)

if(params$ZERO_RESTR){
  ha[ix_params_sys, ix_params_noise] = 0
  ha[ix_params_noise, ix_params_sys] = 0
}

ha_inv = ginv(ha)

Calculate Hessian via optim

hopt = out$hessian

if(params$ZERO_RESTR){
  hopt[ix_params_sys, ix_params_noise] = 0
  hopt[ix_params_noise, ix_params_sys] = 0
}

hopt_inv = ginv(hopt)

3.2 Plot Parameters

Hessian (inverse of information matrix under correct specification) and its normalised version.

In the first column are the

• OPG
• Analytic Hessian
• Hessian obtain from optim

and in the second column their normalised versions.

par(mfcol = c(3,2))

opg %>% corrplot(is.corr = FALSE)
ha %>% corrplot(is.corr = FALSE)
hopt %>% corrplot(is.corr = FALSE)
opg %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt %>% cov2cor() %>% corrplot() %>% suppressWarnings()

3.2.1 System Paramters

par(mfcol = c(3,2))

opg[ix_params_sys, ix_params_sys] %>% corrplot(is.corr = FALSE)
ha[ix_params_sys, ix_params_sys] %>% corrplot(is.corr = FALSE)
hopt[ix_params_sys, ix_params_sys] %>% corrplot(is.corr = FALSE)
opg[ix_params_sys, ix_params_sys] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha[ix_params_sys, ix_params_sys] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt[ix_params_sys, ix_params_sys] %>%cov2cor() %>% corrplot() %>% suppressWarnings()

3.2.2 Noise Parameters

3.2.2.1 Noise Parameters with Strict Zeros

opg[ix_params_sys, ix_params_noise] = 0
opg[ix_params_noise, ix_params_sys] = 0

opg[ix_params_noise_b, ix_params_noise_sgt] = 0
opg[ix_params_noise_sgt, ix_params_noise_b] = 0

opg[ix_params_noise_sgt[1:3], ix_params_noise_sgt[4:6]] = 0
opg[ix_params_noise_sgt[4:6], ix_params_noise_sgt[1:3]] = 0

opg %>% zapsmall()
##           [,1]      [,2]     [,3]     [,4]     [,5]      [,6]     [,7]     [,8]
##  [1,]  8.30260   5.51231 -0.37724  2.80218 -8.21639  -4.26124  2.53252  3.27255
##  [2,]  5.51231  29.31015 -3.27003 -1.30250 -6.63013 -27.34849  1.12524  7.23095
##  [3,] -0.37724  -3.27003  4.37678 -8.81453  1.63817  -0.32804 -0.76047 -0.82736
##  [4,]  2.80218  -1.30250 -8.81453 52.56240 -5.27292  19.03739  2.65373 -2.56562
##  [5,] -8.21639  -6.63013  1.63817 -5.27292  8.75659   4.66961 -2.52694 -3.42208
##  [6,] -4.26124 -27.34849 -0.32804 19.03739  4.66961  32.68539  0.36090 -7.17842
##  [7,]  2.53252   1.12524 -0.76047  2.65373 -2.52694   0.36090  1.27433  1.17358
##  [8,]  3.27255   7.23095 -0.82736 -2.56562 -3.42208  -7.17842  1.17358  3.77695
##  [9,]  0.02605   0.19201 -0.54652  0.30418 -0.15658  -0.05488  0.18029  0.27542
## [10,] -1.47012   1.80855 -0.81566  0.50156  1.07335  -1.91806 -0.52383  0.11658
## [11,]  0.20368   0.02012  0.39924 -0.99184 -0.07523  -0.42292 -0.07180  0.01449
## [12,] -0.48979  -1.34146  0.94565 -3.17236  0.80975  -0.11514 -0.42854 -0.45038
## [13,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [14,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [15,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [16,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [17,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [18,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [19,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [20,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [21,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
## [22,]  0.00000   0.00000  0.00000  0.00000  0.00000   0.00000  0.00000  0.00000
##           [,9]    [,10]    [,11]    [,12]    [,13]    [,14]   [,15]    [,16]
##  [1,]  0.02605 -1.47012  0.20368 -0.48979  0.00000  0.00000 0.00000  0.00000
##  [2,]  0.19201  1.80855  0.02012 -1.34146  0.00000  0.00000 0.00000  0.00000
##  [3,] -0.54652 -0.81566  0.39924  0.94565  0.00000  0.00000 0.00000  0.00000
##  [4,]  0.30418  0.50156 -0.99184 -3.17236  0.00000  0.00000 0.00000  0.00000
##  [5,] -0.15658  1.07335 -0.07523  0.80975  0.00000  0.00000 0.00000  0.00000
##  [6,] -0.05488 -1.91806 -0.42292 -0.11514  0.00000  0.00000 0.00000  0.00000
##  [7,]  0.18029 -0.52383 -0.07180 -0.42854  0.00000  0.00000 0.00000  0.00000
##  [8,]  0.27542  0.11658  0.01449 -0.45038  0.00000  0.00000 0.00000  0.00000
##  [9,]  2.44750  4.33006 -0.44009 -0.53731  0.00000  0.00000 0.00000  0.00000
## [10,]  4.33006 19.07369 -0.67460 -1.99645  0.00000  0.00000 0.00000  0.00000
## [11,] -0.44009 -0.67460  0.30339  0.38985  0.00000  0.00000 0.00000  0.00000
## [12,] -0.53731 -1.99645  0.38985  2.05647  0.00000  0.00000 0.00000  0.00000
## [13,]  0.00000  0.00000  0.00000  0.00000  2.76124  3.10288 1.01643 -0.62951
## [14,]  0.00000  0.00000  0.00000  0.00000  3.10288 26.07433 3.52805 -7.89943
## [15,]  0.00000  0.00000  0.00000  0.00000  1.01643  3.52805 2.29553  3.91247
## [16,]  0.00000  0.00000  0.00000  0.00000 -0.62951 -7.89943 3.91247 31.59374
## [17,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000  0.00000
## [18,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000  0.00000
## [19,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000  0.00000
## [20,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000  0.00000
## [21,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000  0.00000
## [22,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000  0.00000
##         [,17]   [,18] [,19]    [,20]    [,21]    [,22]
##  [1,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
##  [2,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
##  [3,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
##  [4,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
##  [5,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
##  [6,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
##  [7,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
##  [8,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
##  [9,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [10,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [11,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [12,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [13,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [14,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [15,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [16,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [17,] 0.87789 0.07510     0  0.00000  0.00000  0.00000
## [18,] 0.07510 0.13163     0  0.00000  0.00000  0.00000
## [19,] 0.00000 0.00000     0  0.00000  0.00000  0.00000
## [20,] 0.00000 0.00000     0  0.50289 -0.00350 -0.00029
## [21,] 0.00000 0.00000     0 -0.00350  0.05719  0.00166
## [22,] 0.00000 0.00000     0 -0.00029  0.00166  0.00005
opg %>% ginv() %>% zapsmall()
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
##  [1,]  7.74 -0.78 -2.40  0.24  6.79 -0.70 -3.40  0.28 -0.20 -0.02 -1.11 -0.42
##  [2,] -0.78  2.43  0.46 -0.82 -1.73  2.68 -2.64 -0.05  0.09  0.04 -0.46  0.34
##  [3,] -2.40  0.46  1.22 -0.06 -2.16  0.42  0.83  0.04  0.10  0.01 -0.01  0.17
##  [4,]  0.24 -0.82 -0.06  0.33  0.61 -0.93  0.91  0.07  0.00 -0.02  0.18 -0.08
##  [5,]  6.79 -1.73 -2.16  0.61  6.84 -1.86 -0.88  0.25 -0.16 -0.05 -0.41 -0.51
##  [6,] -0.70  2.68  0.42 -0.93 -1.86  3.07 -3.39  0.17  0.07  0.06 -0.55  0.36
##  [7,] -3.40 -2.64  0.83  0.91 -0.88 -3.39  8.91 -1.22 -0.16  0.04  1.49 -0.04
##  [8,]  0.28 -0.05  0.04  0.07  0.25  0.17 -1.22  1.12 -0.02 -0.01 -0.16  0.04
##  [9,] -0.20  0.09  0.10  0.00 -0.16  0.07 -0.16 -0.02  0.97 -0.20  1.21 -0.18
## [10,] -0.02  0.04  0.01 -0.02 -0.05  0.06  0.04 -0.01 -0.20  0.11 -0.18  0.10
## [11,] -1.11 -0.46 -0.01  0.18 -0.41 -0.55  1.49 -0.16  1.21 -0.18  6.85 -1.02
## [12,] -0.42  0.34  0.17 -0.08 -0.51  0.36 -0.04  0.04 -0.18  0.10 -1.02  0.87
## [13,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [14,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [15,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [16,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [17,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [18,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [19,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [20,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [21,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [22,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
##       [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]    [,21]     [,22]
##  [1,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
##  [2,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
##  [3,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
##  [4,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
##  [5,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
##  [6,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
##  [7,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
##  [8,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
##  [9,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
## [10,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
## [11,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
## [12,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
## [13,]  0.48 -0.01 -0.27  0.04  0.00  0.00     0  0.00     0.00      0.00
## [14,] -0.01  0.08 -0.19  0.04  0.00  0.00     0  0.00     0.00      0.00
## [15,] -0.27 -0.19  1.18 -0.20  0.00  0.00     0  0.00     0.00      0.00
## [16,]  0.04  0.04 -0.20  0.07  0.00  0.00     0  0.00     0.00      0.00
## [17,]  0.00  0.00  0.00  0.00  1.20 -0.68     0  0.00     0.00      0.00
## [18,]  0.00  0.00  0.00  0.00 -0.68  7.99     0  0.00     0.00      0.00
## [19,]  0.00  0.00  0.00  0.00  0.00  0.00     0  0.00     0.00      0.00
## [20,]  0.00  0.00  0.00  0.00  0.00  0.00     0  2.01    -1.72     63.53
## [21,]  0.00  0.00  0.00  0.00  0.00  0.00     0 -1.72   158.15  -4856.51
## [22,]  0.00  0.00  0.00  0.00  0.00  0.00     0 63.53 -4856.51 167688.11

hopt[ix_params_sys, ix_params_noise] = 0
hopt[ix_params_noise, ix_params_sys] = 0

hopt[ix_params_noise_b, ix_params_noise_sgt] = 0
hopt[ix_params_noise_sgt, ix_params_noise_b] = 0

hopt[ix_params_noise_sgt[1:3], ix_params_noise_sgt[4:6]] = 0
hopt[ix_params_noise_sgt[4:6], ix_params_noise_sgt[1:3]] = 0

hopt %>% zapsmall()
##           [,1]      [,2]      [,3]      [,4]     [,5]      [,6]     [,7]
##  [1,]  6.80671   1.15824   1.58211  -2.24422 -6.18947  -2.71321  1.22165
##  [2,]  1.15824  24.23756  -4.78356   7.13560 -3.35882 -19.07695  0.95329
##  [3,]  1.58211  -4.78356   6.88489 -14.63398  0.66565  -0.92152 -0.36172
##  [4,] -2.24422   7.13560 -14.63398  68.70189 -2.18633  15.49050  0.12931
##  [5,] -6.18947  -3.35882   0.66565  -2.18633  6.86046   3.24671 -1.21063
##  [6,] -2.71321 -19.07695  -0.92152  15.49050  3.24671  23.47750  0.20038
##  [7,]  1.22165   0.95329  -0.36172   0.12931 -1.21063   0.20038  1.20988
##  [8,]  1.57283   3.63632  -0.77461  -2.12066 -1.63627  -3.87901  0.60253
##  [9,] -1.80201  -1.05038  -1.76138   3.38608  1.54798   2.45213 -0.21297
## [10,]  0.04229   0.60600  -0.51997  -1.46983  1.00414   0.14682  1.22563
## [11,]  0.51867  -0.39122   0.44157  -1.27603 -0.43692  -0.06173  0.10368
## [12,] -0.31442   0.44981   0.12153  -1.27525  0.00428  -1.01057 -0.29194
## [13,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [14,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [15,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [16,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [17,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [18,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [19,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [20,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [21,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [22,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
##           [,8]     [,9]    [,10]    [,11]    [,12]    [,13]    [,14]   [,15]
##  [1,]  1.57283 -1.80201  0.04229  0.51867 -0.31442  0.00000  0.00000 0.00000
##  [2,]  3.63632 -1.05038  0.60600 -0.39122  0.44981  0.00000  0.00000 0.00000
##  [3,] -0.77461 -1.76138 -0.51997  0.44157  0.12153  0.00000  0.00000 0.00000
##  [4,] -2.12066  3.38608 -1.46983 -1.27603 -1.27525  0.00000  0.00000 0.00000
##  [5,] -1.63627  1.54798  1.00414 -0.43692  0.00428  0.00000  0.00000 0.00000
##  [6,] -3.87901  2.45213  0.14682 -0.06173 -1.01057  0.00000  0.00000 0.00000
##  [7,]  0.60253 -0.21297  1.22563  0.10368 -0.29194  0.00000  0.00000 0.00000
##  [8,]  2.54610  0.53657  0.80833  0.05397  0.61825  0.00000  0.00000 0.00000
##  [9,]  0.53657  4.84800  4.80466 -0.90319 -0.64090  0.00000  0.00000 0.00000
## [10,]  0.80833  4.80466 21.27074 -0.58546 -3.19909  0.00000  0.00000 0.00000
## [11,]  0.05397 -0.90319 -0.58546  0.36647  0.29328  0.00000  0.00000 0.00000
## [12,]  0.61825 -0.64090 -3.19909  0.29328  1.87076  0.00000  0.00000 0.00000
## [13,]  0.00000  0.00000  0.00000  0.00000  0.00000  2.65480  0.59093 0.66053
## [14,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.59093 19.30288 3.20398
## [15,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.66053  3.20398 3.75069
## [16,]  0.00000  0.00000  0.00000  0.00000  0.00000 -1.06378 -5.10908 2.86177
## [17,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [18,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [19,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [20,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [21,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [22,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
##          [,16]    [,17]    [,18] [,19]    [,20]    [,21]    [,22]
##  [1,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [2,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [3,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [4,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [5,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [6,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [7,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [8,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [9,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [10,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [11,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [12,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [13,] -1.06378  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [14,] -5.10908  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [15,]  2.86177  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [16,] 31.95176  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [17,]  0.00000  2.46573 -0.14302     0  0.00000  0.00000  0.00000
## [18,]  0.00000 -0.14302  0.17574     0  0.00000  0.00000  0.00000
## [19,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [20,]  0.00000  0.00000  0.00000     0  0.44235 -0.05234 -0.00107
## [21,]  0.00000  0.00000  0.00000     0 -0.05234  0.07686  0.00238
## [22,]  0.00000  0.00000  0.00000     0 -0.00107  0.00238  0.00008
hopt %>% ginv() %>% zapsmall()
##        [,1]   [,2]  [,3]  [,4]  [,5]   [,6]   [,7]  [,8]  [,9] [,10] [,11]
##  [1,]  4.03   0.83 -1.47 -0.36  3.18   0.64  -0.93 -1.50  0.44 -0.08  0.25
##  [2,]  0.83   7.33  0.06 -2.65 -2.12   8.70 -11.30  1.80 -2.22  0.58 -5.32
##  [3,] -1.47   0.06  0.91  0.09 -1.25   0.15  -0.08  0.74 -0.14  0.03 -0.19
##  [4,] -0.36  -2.65  0.09  1.00  0.75  -3.18   4.13 -0.61  0.84 -0.22  2.09
##  [5,]  3.18  -2.12 -1.25  0.75  3.84  -2.86   3.76 -1.86  1.25 -0.33  2.67
##  [6,]  0.64   8.70  0.15 -3.18 -2.86  10.49 -13.66  2.50 -2.91  0.74 -6.93
##  [7,] -0.93 -11.30 -0.08  4.13  3.76 -13.66  19.06 -3.53  3.94 -1.04  8.86
##  [8,] -1.50   1.80  0.74 -0.61 -1.86   2.50  -3.53  2.37 -1.33  0.22 -2.37
##  [9,]  0.44  -2.22 -0.14  0.84  1.25  -2.91   3.94 -1.33  1.75 -0.39  3.89
## [10,] -0.08   0.58  0.03 -0.22 -0.33   0.74  -1.04  0.22 -0.39  0.16 -0.89
## [11,]  0.25  -5.32 -0.19  2.09  2.67  -6.93   8.86 -2.37  3.89 -0.89 13.07
## [12,]  1.00  -0.02 -0.40 -0.03  0.73  -0.08   0.33 -0.84  0.00  0.15 -1.05
## [13,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [14,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [15,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [16,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [17,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [18,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [19,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [20,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [21,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
## [22,]  0.00   0.00  0.00  0.00  0.00   0.00   0.00  0.00  0.00  0.00  0.00
##       [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]   [,20]    [,21]
##  [1,]  1.00  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
##  [2,] -0.02  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
##  [3,] -0.40  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
##  [4,] -0.03  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
##  [5,]  0.73  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
##  [6,] -0.08  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
##  [7,]  0.33  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
##  [8,] -0.84  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
##  [9,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
## [10,]  0.15  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
## [11,] -1.05  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
## [12,]  1.42  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
## [13,]  0.00  0.41  0.01 -0.10  0.02  0.00  0.00     0    0.00     0.00
## [14,]  0.00  0.01  0.07 -0.07  0.02  0.00  0.00     0    0.00     0.00
## [15,]  0.00 -0.10 -0.07  0.39 -0.05  0.00  0.00     0    0.00     0.00
## [16,]  0.00  0.02  0.02 -0.05  0.04  0.00  0.00     0    0.00     0.00
## [17,]  0.00  0.00  0.00  0.00  0.00  0.43  0.35     0    0.00     0.00
## [18,]  0.00  0.00  0.00  0.00  0.00  0.35  5.97     0    0.00     0.00
## [19,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00
## [20,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00     0    2.86    11.08
## [21,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00     0   11.08   236.46
## [22,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00     0 -294.85 -6967.65
##           [,22]
##  [1,]      0.00
##  [2,]      0.00
##  [3,]      0.00
##  [4,]      0.00
##  [5,]      0.00
##  [6,]      0.00
##  [7,]      0.00
##  [8,]      0.00
##  [9,]      0.00
## [10,]      0.00
## [11,]      0.00
## [12,]      0.00
## [13,]      0.00
## [14,]      0.00
## [15,]      0.00
## [16,]      0.00
## [17,]      0.00
## [18,]      0.00
## [19,]      0.00
## [20,]   -294.85
## [21,]  -6967.65
## [22,] 218382.47

ha[ix_params_sys, ix_params_noise] = 0
ha[ix_params_noise, ix_params_sys] = 0

ha[ix_params_noise_b, ix_params_noise_sgt] = 0
ha[ix_params_noise_sgt, ix_params_noise_b] = 0

ha[ix_params_noise_sgt[1:3], ix_params_noise_sgt[4:6]] = 0
ha[ix_params_noise_sgt[4:6], ix_params_noise_sgt[1:3]] = 0

ha %>% zapsmall()
##           [,1]      [,2]      [,3]      [,4]     [,5]      [,6]     [,7]
##  [1,]  6.64779   1.38285   1.55838   0.21292 -6.13935  -2.15539  1.37702
##  [2,]  1.38285  23.29748  -4.65762   4.68690 -2.95251 -18.50730  0.22495
##  [3,]  1.55838  -4.65762   7.16766 -11.52276  0.54254  -0.78910 -0.36629
##  [4,]  0.21292   4.68690 -11.52276  62.16895 -3.75032  14.43507  0.44700
##  [5,] -6.13935  -2.95251   0.54254  -3.75032  6.84693   2.99103 -1.12235
##  [6,] -2.15539 -18.50730  -0.78910  14.43507  2.99103  23.52999  0.63428
##  [7,]  1.37702   0.22495  -0.36629   0.44700 -1.12235   0.63428  1.04505
##  [8,]  1.57400   3.65093  -0.85760  -3.95510 -1.63601  -4.14276  0.63166
##  [9,] -1.29375  -0.69725  -1.13021  -0.37457  1.34147   1.16890 -0.09826
## [10,]  0.14462   0.95949  -0.07914  -1.73853  1.06119  -0.05290  1.15974
## [11,]  0.34404  -0.55196   0.23172  -0.65132 -0.34884   0.14036  0.05165
## [12,] -0.31968   0.50839   0.02581  -1.50673  0.00026  -0.95949 -0.27621
## [13,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [14,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [15,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [16,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [17,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [18,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [19,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [20,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [21,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
## [22,]  0.00000   0.00000   0.00000   0.00000  0.00000   0.00000  0.00000
##           [,8]     [,9]    [,10]    [,11]    [,12]    [,13]    [,14]   [,15]
##  [1,]  1.57400 -1.29375  0.14462  0.34404 -0.31968  0.00000  0.00000 0.00000
##  [2,]  3.65093 -0.69725  0.95949 -0.55196  0.50839  0.00000  0.00000 0.00000
##  [3,] -0.85760 -1.13021 -0.07914  0.23172  0.02581  0.00000  0.00000 0.00000
##  [4,] -3.95510 -0.37457 -1.73853 -0.65132 -1.50673  0.00000  0.00000 0.00000
##  [5,] -1.63601  1.34147  1.06119 -0.34884  0.00026  0.00000  0.00000 0.00000
##  [6,] -4.14276  1.16890 -0.05290  0.14036 -0.95949  0.00000  0.00000 0.00000
##  [7,]  0.63166 -0.09826  1.15974  0.05165 -0.27621  0.00000  0.00000 0.00000
##  [8,]  2.53165  0.54167  0.71112  0.05541  0.61261  0.00000  0.00000 0.00000
##  [9,]  0.54167  4.18309  4.79907 -0.75394 -0.64415  0.00000  0.00000 0.00000
## [10,]  0.71112  4.79907 21.11990 -0.57324 -3.14332  0.00000  0.00000 0.00000
## [11,]  0.05541 -0.75394 -0.57324  0.34160  0.29750  0.00000  0.00000 0.00000
## [12,]  0.61261 -0.64415 -3.14332  0.29750  1.87748  0.00000  0.00000 0.00000
## [13,]  0.00000  0.00000  0.00000  0.00000  0.00000  2.61700  0.71211 0.61947
## [14,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.71211 19.33209 2.95817
## [15,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.61947  2.95817 3.69067
## [16,]  0.00000  0.00000  0.00000  0.00000  0.00000 -1.24011 -5.22227 2.89039
## [17,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [18,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [19,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [20,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [21,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
## [22,]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000 0.00000
##          [,16]    [,17]    [,18] [,19]    [,20]    [,21]    [,22]
##  [1,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [2,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [3,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [4,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [5,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [6,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [7,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [8,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
##  [9,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [10,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [11,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [12,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [13,] -1.24011  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [14,] -5.22227  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [15,]  2.89039  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [16,] 32.05416  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [17,]  0.00000  2.20306 -0.10606     0  0.00000  0.00000  0.00000
## [18,]  0.00000 -0.10606  0.18180     0  0.00000  0.00000  0.00000
## [19,]  0.00000  0.00000  0.00000     0  0.00000  0.00000  0.00000
## [20,]  0.00000  0.00000  0.00000     0  0.44234 -0.05233 -0.00107
## [21,]  0.00000  0.00000  0.00000     0 -0.05233  0.07686  0.00238
## [22,]  0.00000  0.00000  0.00000     0 -0.00107  0.00238  0.00008
ha %>% ginv() %>% zapsmall()
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11] [,12]
##  [1,]  6.90 -1.24 -3.02  0.13  5.92 -1.57 -0.83 -2.27 -0.05  0.07 -0.97  1.72
##  [2,] -1.24  0.72  0.69 -0.15 -1.22  0.75 -0.23  0.29  0.18 -0.05  0.77 -0.43
##  [3,] -3.02  0.69  1.67 -0.01 -2.54  0.81  0.31  1.21  0.20 -0.09  0.90 -0.89
##  [4,]  0.13 -0.15 -0.01  0.08  0.21 -0.18  0.12  0.08  0.01 -0.01  0.05  0.02
##  [5,]  5.92 -1.22 -2.54  0.21  5.44 -1.56 -0.35 -1.70 -0.05  0.02 -0.56  1.36
##  [6,] -1.57  0.75  0.81 -0.18 -1.56  0.94 -0.53  0.69 -0.06 -0.01  0.21 -0.51
##  [7,] -0.83 -0.23  0.31  0.12 -0.35 -0.53  2.86 -0.82  0.64 -0.17  1.00  0.20
##  [8,] -2.27  0.29  1.21  0.08 -1.70  0.69 -0.82  2.57 -0.62  0.02 -0.59 -1.11
##  [9,] -0.05  0.18  0.20  0.01 -0.05 -0.06  0.64 -0.62  1.11 -0.23  2.43 -0.17
## [10,]  0.07 -0.05 -0.09 -0.01  0.02 -0.01 -0.17  0.02 -0.23  0.12 -0.54  0.20
## [11,] -0.97  0.77  0.90  0.05 -0.56  0.21  1.00 -0.59  2.43 -0.54  9.68 -1.50
## [12,]  1.72 -0.43 -0.89  0.02  1.36 -0.51  0.20 -1.11 -0.17  0.20 -1.50  1.62
## [13,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [14,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [15,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [16,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [17,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [18,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [19,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [20,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [21,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
## [22,]  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00
##       [,13] [,14] [,15] [,16] [,17] [,18] [,19]   [,20]    [,21]     [,22]
##  [1,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
##  [2,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
##  [3,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
##  [4,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
##  [5,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
##  [6,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
##  [7,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
##  [8,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
##  [9,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
## [10,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
## [11,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
## [12,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
## [13,]  0.41  0.01 -0.09  0.03  0.00  0.00     0    0.00     0.00      0.00
## [14,]  0.01  0.07 -0.07  0.02  0.00  0.00     0    0.00     0.00      0.00
## [15,] -0.09 -0.07  0.38 -0.05  0.00  0.00     0    0.00     0.00      0.00
## [16,]  0.03  0.02 -0.05  0.04  0.00  0.00     0    0.00     0.00      0.00
## [17,]  0.00  0.00  0.00  0.00  0.47  0.27     0    0.00     0.00      0.00
## [18,]  0.00  0.00  0.00  0.00  0.27  5.66     0    0.00     0.00      0.00
## [19,]  0.00  0.00  0.00  0.00  0.00  0.00     0    0.00     0.00      0.00
## [20,]  0.00  0.00  0.00  0.00  0.00  0.00     0    2.86    11.08   -294.85
## [21,]  0.00  0.00  0.00  0.00  0.00  0.00     0   11.08   236.46  -6967.70
## [22,]  0.00  0.00  0.00  0.00  0.00  0.00     0 -294.85 -6967.70 218384.10

c(opg[ix_params_noise, ix_params_noise], 
  ha[ix_params_noise, ix_params_noise], 
  hopt[ix_params_noise, ix_params_noise]) %>% 
  array(dim = c(10, 10, 3)) %>% 
  zapsmall()
## , , 1
## 
##           [,1]     [,2]    [,3]     [,4]    [,5]    [,6] [,7]     [,8]     [,9]
##  [1,]  2.76124  3.10288 1.01643 -0.62951 0.00000 0.00000    0  0.00000  0.00000
##  [2,]  3.10288 26.07433 3.52805 -7.89943 0.00000 0.00000    0  0.00000  0.00000
##  [3,]  1.01643  3.52805 2.29553  3.91247 0.00000 0.00000    0  0.00000  0.00000
##  [4,] -0.62951 -7.89943 3.91247 31.59374 0.00000 0.00000    0  0.00000  0.00000
##  [5,]  0.00000  0.00000 0.00000  0.00000 0.87789 0.07510    0  0.00000  0.00000
##  [6,]  0.00000  0.00000 0.00000  0.00000 0.07510 0.13163    0  0.00000  0.00000
##  [7,]  0.00000  0.00000 0.00000  0.00000 0.00000 0.00000    0  0.00000  0.00000
##  [8,]  0.00000  0.00000 0.00000  0.00000 0.00000 0.00000    0  0.50289 -0.00350
##  [9,]  0.00000  0.00000 0.00000  0.00000 0.00000 0.00000    0 -0.00350  0.05719
## [10,]  0.00000  0.00000 0.00000  0.00000 0.00000 0.00000    0 -0.00029  0.00166
##          [,10]
##  [1,]  0.00000
##  [2,]  0.00000
##  [3,]  0.00000
##  [4,]  0.00000
##  [5,]  0.00000
##  [6,]  0.00000
##  [7,]  0.00000
##  [8,] -0.00029
##  [9,]  0.00166
## [10,]  0.00005
## 
## , , 2
## 
##           [,1]     [,2]    [,3]     [,4]     [,5]     [,6] [,7]     [,8]
##  [1,]  2.61700  0.71211 0.61947 -1.24011  0.00000  0.00000    0  0.00000
##  [2,]  0.71211 19.33209 2.95817 -5.22227  0.00000  0.00000    0  0.00000
##  [3,]  0.61947  2.95817 3.69067  2.89039  0.00000  0.00000    0  0.00000
##  [4,] -1.24011 -5.22227 2.89039 32.05416  0.00000  0.00000    0  0.00000
##  [5,]  0.00000  0.00000 0.00000  0.00000  2.20306 -0.10606    0  0.00000
##  [6,]  0.00000  0.00000 0.00000  0.00000 -0.10606  0.18180    0  0.00000
##  [7,]  0.00000  0.00000 0.00000  0.00000  0.00000  0.00000    0  0.00000
##  [8,]  0.00000  0.00000 0.00000  0.00000  0.00000  0.00000    0  0.44234
##  [9,]  0.00000  0.00000 0.00000  0.00000  0.00000  0.00000    0 -0.05233
## [10,]  0.00000  0.00000 0.00000  0.00000  0.00000  0.00000    0 -0.00107
##           [,9]    [,10]
##  [1,]  0.00000  0.00000
##  [2,]  0.00000  0.00000
##  [3,]  0.00000  0.00000
##  [4,]  0.00000  0.00000
##  [5,]  0.00000  0.00000
##  [6,]  0.00000  0.00000
##  [7,]  0.00000  0.00000
##  [8,] -0.05233 -0.00107
##  [9,]  0.07686  0.00238
## [10,]  0.00238  0.00008
## 
## , , 3
## 
##           [,1]     [,2]    [,3]     [,4]     [,5]     [,6] [,7]     [,8]
##  [1,]  2.65480  0.59093 0.66053 -1.06378  0.00000  0.00000    0  0.00000
##  [2,]  0.59093 19.30288 3.20398 -5.10908  0.00000  0.00000    0  0.00000
##  [3,]  0.66053  3.20398 3.75069  2.86177  0.00000  0.00000    0  0.00000
##  [4,] -1.06378 -5.10908 2.86177 31.95176  0.00000  0.00000    0  0.00000
##  [5,]  0.00000  0.00000 0.00000  0.00000  2.46573 -0.14302    0  0.00000
##  [6,]  0.00000  0.00000 0.00000  0.00000 -0.14302  0.17574    0  0.00000
##  [7,]  0.00000  0.00000 0.00000  0.00000  0.00000  0.00000    0  0.00000
##  [8,]  0.00000  0.00000 0.00000  0.00000  0.00000  0.00000    0  0.44235
##  [9,]  0.00000  0.00000 0.00000  0.00000  0.00000  0.00000    0 -0.05234
## [10,]  0.00000  0.00000 0.00000  0.00000  0.00000  0.00000    0 -0.00107
##           [,9]    [,10]
##  [1,]  0.00000  0.00000
##  [2,]  0.00000  0.00000
##  [3,]  0.00000  0.00000
##  [4,]  0.00000  0.00000
##  [5,]  0.00000  0.00000
##  [6,]  0.00000  0.00000
##  [7,]  0.00000  0.00000
##  [8,] -0.05234 -0.00107
##  [9,]  0.07686  0.00238
## [10,]  0.00238  0.00008

hopt[ix_params_noise, ix_params_noise] %>% ginv() %>% zapsmall()
##        [,1]  [,2]  [,3]  [,4] [,5] [,6] [,7]    [,8]     [,9]     [,10]
##  [1,]  0.41  0.01 -0.10  0.02 0.00 0.00    0    0.00     0.00      0.00
##  [2,]  0.01  0.07 -0.07  0.02 0.00 0.00    0    0.00     0.00      0.00
##  [3,] -0.10 -0.07  0.39 -0.05 0.00 0.00    0    0.00     0.00      0.00
##  [4,]  0.02  0.02 -0.05  0.04 0.00 0.00    0    0.00     0.00      0.00
##  [5,]  0.00  0.00  0.00  0.00 0.43 0.35    0    0.00     0.00      0.00
##  [6,]  0.00  0.00  0.00  0.00 0.35 5.97    0    0.00     0.00      0.00
##  [7,]  0.00  0.00  0.00  0.00 0.00 0.00    0    0.00     0.00      0.00
##  [8,]  0.00  0.00  0.00  0.00 0.00 0.00    0    2.86    11.08   -294.85
##  [9,]  0.00  0.00  0.00  0.00 0.00 0.00    0   11.08   236.46  -6967.65
## [10,]  0.00  0.00  0.00  0.00 0.00 0.00    0 -294.85 -6967.65 218382.47

hopt[ix_params_noise_sgt, ix_params_noise_sgt] %>% ginv() %>% zapsmall()
##      [,1] [,2]       [,3]   [,4]    [,5]     [,6]
## [1,]  0.4  0.3        0.4    0.0     0.0      0.0
## [2,]  0.3  6.0        3.4    0.0     0.0      0.0
## [3,]  0.4  3.4 -2180917.9    0.0     0.0      0.0
## [4,]  0.0  0.0        0.0    2.9    11.1   -294.9
## [5,]  0.0  0.0        0.0   11.1   236.5  -6967.7
## [6,]  0.0  0.0        0.0 -294.9 -6967.7 218382.5

eigen(hopt[ix_params_noise, ix_params_noise])$values
##  [1]  3.389182e+01  1.856179e+01  3.248406e+00  2.474626e+00  1.958111e+00
##  [6]  4.496991e-01  1.668400e-01  6.958517e-02  4.574457e-06 -4.585225e-07

eigen(hopt[ix_params_noise, ix_params_noise])$values
##  [1]  3.389182e+01  1.856179e+01  3.248406e+00  2.474626e+00  1.958111e+00
##  [6]  4.496991e-01  1.668400e-01  6.958517e-02  4.574457e-06 -4.585225e-07

eigen(hopt[ix_params_noise_b, ix_params_noise_b])$values
## [1] 33.891818 18.561787  3.248406  1.958111

eigen(hopt[ix_params_noise_sgt, ix_params_noise_sgt])$values
## [1]  2.474626e+00  4.496991e-01  1.668400e-01  6.958517e-02  4.574457e-06
## [6] -4.585225e-07

eigen(hopt[ix_params_noise_sgt[1:3], ix_params_noise_sgt[1:3]])$values
## [1]  2.474626e+00  1.668400e-01 -4.585225e-07

hopt[ix_params_noise_sgt[1:3], ix_params_noise_sgt[1:3]] %>% zapsmall()
##            [,1]       [,2]   [,3]
## [1,]  2.4657285 -0.1430155  3e-07
## [2,] -0.1430155  0.1757371  2e-07
## [3,]  0.0000003  0.0000002 -5e-07
hopt[ix_params_noise_sgt[1:3], ix_params_noise_sgt[1:3]] %>% ginv() %>% zapsmall()
##      [,1] [,2]       [,3]
## [1,]  0.4  0.3        0.4
## [2,]  0.3  6.0        3.4
## [3,]  0.4  3.4 -2180917.9

eigen(hopt[ix_params_noise_sgt[4:6], ix_params_noise_sgt[4:6]])$values
## [1] 4.496991e-01 6.958517e-02 4.574457e-06

hopt[ix_params_noise_sgt[4:6], ix_params_noise_sgt[4:6]] %>% zapsmall()
##            [,1]       [,2]       [,3]
## [1,]  0.4423478 -0.0523373 -0.0010726
## [2,] -0.0523373  0.0768619  0.0023817
## [3,] -0.0010726  0.0023817  0.0000791
hopt[ix_params_noise_sgt[4:6], ix_params_noise_sgt[4:6]] %>% ginv() %>% zapsmall()
##         [,1]     [,2]      [,3]
## [1,]    2.86    11.08   -294.85
## [2,]   11.08   236.46  -6967.65
## [3,] -294.85 -6967.65 218382.47

par(mfcol = c(3,2))

opg[ix_params_noise, ix_params_noise] %>% corrplot(is.corr = FALSE)
ha[ix_params_noise, ix_params_noise] %>% corrplot(is.corr = FALSE)
hopt[ix_params_noise, ix_params_noise] %>% corrplot(is.corr = FALSE)
opg[ix_params_noise, ix_params_noise] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha[ix_params_noise, ix_params_noise] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt[ix_params_noise, ix_params_noise] %>%cov2cor() %>% corrplot() %>% suppressWarnings()

B matrix

par(mfcol = c(3,2))

opg[ix_params_noise_b, ix_params_noise_b] %>% corrplot(is.corr = FALSE)
ha[ix_params_noise_b, ix_params_noise_b] %>% corrplot(is.corr = FALSE)
hopt[ix_params_noise_b, ix_params_noise_b] %>% corrplot(is.corr = FALSE)
opg[ix_params_noise_b, ix_params_noise_b] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha[ix_params_noise_b, ix_params_noise_b] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt[ix_params_noise_b, ix_params_noise_b] %>%cov2cor() %>% corrplot() %>% suppressWarnings()

Distribution parameters

par(mfcol = c(3,2))

opg[ix_params_noise_sgt, ix_params_noise_sgt] %>% corrplot(is.corr = FALSE)
ha[ix_params_noise_sgt, ix_params_noise_sgt] %>% corrplot(is.corr = FALSE)
hopt[ix_params_noise_sgt, ix_params_noise_sgt] %>% corrplot(is.corr = FALSE)
opg[ix_params_noise_sgt, ix_params_noise_sgt] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
ha[ix_params_noise_sgt, ix_params_noise_sgt] %>% cov2cor() %>% corrplot() %>% suppressWarnings()
hopt[ix_params_noise_sgt, ix_params_noise_sgt] %>%cov2cor() %>% corrplot() %>% suppressWarnings()

3.3 Diagonal elemets

We compare the diagonal elements of

• the OPG matrix
• Hessian obtained from optim
• analytic Hessian

3.3.1 Quantitative

The columns starting with diff_ contain the absolute difference between the different objects

(
  comp_diag = tibble(name = c(rep("AR", 4),
                              rep("MA_bwd", 4),
                              rep("MA_fwd", 4),
                              rep("B", 4),
                              rep("distr", 6)),
                     type = c(rep("system", 12), 
                              rep("noise", 10)),
                     a_opg_diag = diag(opg),
                     b_ha_diag = diag(ha),
                     c_hopt_diag = diag(hopt))
)

comp_diag %>% 
  dtable()
## This version of bslib is designed to work with rmarkdown version 2.7 or higher.

We need the following helper function which calculates the sum of the absolute differences between the elements pertaining to (OPG, Hessian optim, Hessian analytic)

obtain_diff_abs_mat = function(df){
  
  mat  = df %>% 
    as.matrix()
  
  names_df = colnames(df)
  
  res = matrix(nrow = ncol(df), ncol = ncol(df))
  
  for (ix_row in 1:ncol(df)){
    for (ix_col in 1:ncol(df)) {
      res[ix_row, ix_col] = (mat[, ix_row] - mat[, ix_col]) %>% abs() %>% sum(na.rm = TRUE)
    }
  }
  
  rownames(res) = names_df
  colnames(res) = names_df
  
  return(res)
}

Investigate the differences between the diagonal elements for each group (AR, MA backwards, MA forwards, B, distribution).

comp_diag %>% select(contains("diag")) %>% obtain_diff_abs_mat()
##             a_opg_diag b_ha_diag c_hopt_diag
## a_opg_diag      0.0000 46.801195   52.865005
## b_ha_diag      46.8012  0.000000    9.505489
## c_hopt_diag    52.8650  9.505489    0.000000

3.3.2 Plots

Plot AR parameters

comp_diag %>% 
  filter(name == "AR") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Diagonal Elements for System Parameters: AR")

Plot MA backwards parameters

comp_diag %>% 
  filter(name == "MA_bwd") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Diagonal Elements for System Parameters: MA backwards")

Plot MA forwards parameters

comp_diag %>% 
  filter(name == "MA_fwd") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Diagonal Elements for System Parameters: MA forwards")

Plot noise parameters

comp_diag %>% 
  filter(name == "B") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Diagonal Elements for Noise Parameters: B")

comp_diag %>% 
  filter(name == "distr") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Diagonal Elements for Noise Parameters: Distribution")

hopt %>% diag() %>% sqrt()
## Warning in sqrt(.): NaNs produced
##  [1] 2.608967522 4.923165175 2.623907143 8.288660201 2.619247677 4.845358745
##  [7] 1.099944854 1.595651871 2.201818164 4.612020757 0.605369915 1.367756995
## [13] 1.629354911 4.393503760 1.936669761 5.652588739 1.570263821 0.419210130
## [19]         NaN 0.665092310 0.277239802 0.008894932

4 Compare Standard Deviations

4.1 Quantitative

Difference between all possible options for calculating the standard deviations of the parameters.

comp_sd = tibble(name = c(rep("AR", 4),
                          rep("MA_bwd", 4),
                          rep("MA_fwd", 4),
                          rep("B", 4),
                          rep("distr", 6)),
                 type = c(rep("system", 12), 
                          rep("noise", 10))) %>% 
  unite("type_name", type, name, remove = FALSE) %>% 
  mutate(aa_opg_inv  = diag(ginv(opg)) %>% sqrt(),
         ab_hopt_inv = diag(ginv(hopt)) %>% sqrt(),
         ac_ha_inv   = diag(ginv(ha)) %>% sqrt(),
         ba_rob_hopt = diag(ginv(hopt) %*% opg %*% ginv(hopt)) %>% sqrt(),
         bb_rob_ha   = diag(ginv(ha) %*% opg %*% ginv(ha)) %>% sqrt())

comp_sd %>% 
  mutate(across(where(is.numeric), ~ round(.x, digits = 3))) %>% 
  dtable()

comp_sd %>% 
  select(starts_with("a"), starts_with("b")) %>% 
  obtain_diff_abs_mat()
##             aa_opg_inv ab_hopt_inv  ac_ha_inv ba_rob_hopt bb_rob_ha
## aa_opg_inv     0.00000   70.489367  67.532362   249.89354 231.34195
## ab_hopt_inv   70.48937    0.000000   9.927695   182.59183 172.93476
## ac_ha_inv     67.53236    9.927695   0.000000   189.43229 166.10751
## ba_rob_hopt  249.89354  182.591831 189.432286     0.00000  31.76549
## bb_rob_ha    231.34195  172.934764 166.107507    31.76549   0.00000

comp_sd %>% 
  filter(type == "system") %>% 
  select(starts_with("a"), starts_with("b")) %>% 
  obtain_diff_abs_mat()
##             aa_opg_inv ab_hopt_inv ac_ha_inv ba_rob_hopt bb_rob_ha
## aa_opg_inv    0.000000    8.175851  5.179494    32.95936  14.29749
## ab_hopt_inv   8.175851    0.000000  9.816698    27.93706  18.16038
## ac_ha_inv     5.179494    9.816698  0.000000    34.80916  11.37410
## ba_rob_hopt  32.959362   27.937064 34.809156     0.00000  31.57388
## bb_rob_ha    14.297493   18.160381 11.374097    31.57388   0.00000

comp_sd %>% 
  filter(type == "noise") %>% 
  select(starts_with("a"), starts_with("b")) %>% 
  obtain_diff_abs_mat()
##             aa_opg_inv ab_hopt_inv   ac_ha_inv ba_rob_hopt   bb_rob_ha
## aa_opg_inv     0.00000  62.3135158  62.3528677 216.9341786 217.0444585
## ab_hopt_inv   62.31352   0.0000000   0.1109974 154.6547665 154.7743831
## ac_ha_inv     62.35287   0.1109974   0.0000000 154.6231304 154.7334103
## ba_rob_hopt  216.93418 154.6547665 154.6231304   0.0000000   0.1916113
## bb_rob_ha    217.04446 154.7743831 154.7334103   0.1916113   0.0000000

comp_sd %>% 
  filter(name == "B") %>% 
  select(starts_with("a"), starts_with("b")) %>% 
  obtain_diff_abs_mat()
##             aa_opg_inv ab_hopt_inv  ac_ha_inv ba_rob_hopt  bb_rob_ha
## aa_opg_inv   0.0000000  0.59969563 0.60340252  0.85979233 0.85665695
## ab_hopt_inv  0.5996956  0.00000000 0.01339046  0.29419986 0.30039404
## ac_ha_inv    0.6034025  0.01339046 0.00000000  0.29820929 0.29507391
## ba_rob_hopt  0.8597923  0.29419986 0.29820929  0.00000000 0.02700318
## bb_rob_ha    0.8566569  0.30039404 0.29507391  0.02700318 0.00000000

comp_sd %>% 
  filter(name == "distr") %>% 
  select(starts_with("a"), starts_with("b")) %>% 
  obtain_diff_abs_mat()
##             aa_opg_inv  ab_hopt_inv    ac_ha_inv ba_rob_hopt   bb_rob_ha
## aa_opg_inv     0.00000  61.71382021  61.74946520 216.0743863 216.1878016
## ab_hopt_inv   61.71382   0.00000000   0.09760692 154.3605666 154.4739890
## ac_ha_inv     61.74947   0.09760692   0.00000000 154.3249211 154.4383364
## ba_rob_hopt  216.07439 154.36056662 154.32492111   0.0000000   0.1646081
## bb_rob_ha    216.18780 154.47398903 154.43833639   0.1646081   0.0000000

4.2 Plots

comp_sd %>% 
  filter(name == "AR") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type, -type_name) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Standard Errors of AR Parameters")

comp_sd %>% 
  filter(name == "MA_bwd") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type, -type_name) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Standard Errors of MA backwards Parameters")

comp_sd %>% 
  filter(name == "MA_fwd") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type, -type_name) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Standard Errors of MA forwards Parameters")

comp_sd %>% 
  filter(name == "B") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type, -type_name) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Standard Errors of B Parameters")

PARAMS[ix_params_noise_sgt]
## [1]    -0.5243178     1.4924937 85616.3114036     0.1492405     1.9211332
## [6]     7.8907217

comp_sd %>% 
  filter(name == "distr") %>% 
  rowid_to_column() %>% 
  select(-contains("diff")) %>% 
  select(-name, -type, -type_name) %>% 
  pivot_longer(-rowid) %>% 
  ggplot(aes(x = rowid, y = value, fill = name)) + 
  geom_bar(stat = "identity", position = position_dodge()) + 
  ggtitle("Standard Errors of Distribution Parameters")

5 Results

We choose the OPG for standard deviations in the correctly specified setting and the robust version calculated with analytic Hessian in the misspecified setting

opt_mod_params = fill_tmpl_whf_rev(theta = PARAMS,
                                   tmpl = TMPL)
opt_mod_sd_correctlyspecified = fill_tmpl_whf_rev(theta = comp_sd %>% pull(aa_opg_inv),
                                                  tmpl = TMPL)
opt_mod_sd_misspecified = fill_tmpl_whf_rev(theta = comp_sd %>% pull(ba_rob_hopt),
                                            tmpl = TMPL)

Print AR parameters with non-robust and robust standard errors in wide matrix format.

opt_mod_params$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]   [,3]   [,4]
## [1,]    1    0 -0.107 -0.277
## [2,]    0    1  0.492 -1.013
opt_mod_sd_correctlyspecified$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 2.782 1.104
## [2,]    0    1 1.559 0.575
opt_mod_sd_misspecified$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 2.921 1.308
## [2,]    0    1 6.327 2.331

Print MA backwards parameters with non-robust and robust standard errors in wide matrix format.

opt_mod_params$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 0.065 0.749
## [2,]    0    1 0.259 0.002
opt_mod_sd_correctlyspecified$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 2.615 2.985
## [2,]    0    1 1.751 1.057
opt_mod_sd_misspecified$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]   [,4]
## [1,]    1    0 3.689 10.528
## [2,]    0    1 7.864  3.656

Print MA forwards parameters with non-robust and robust standard errors in wide matrix format.

opt_mod_params$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]   [,3]   [,4]
## [1,]    1    0  0.561  1.671
## [2,]    0    1 -0.241 -0.296
opt_mod_sd_correctlyspecified$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 0.985 2.617
## [2,]    0    1 0.328 0.932
opt_mod_sd_misspecified$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 3.203 7.564
## [2,]    0    1 0.780 2.077

Print static matrix B with non-robust and robust standard errors.

opt_mod_params$B %>% matrix(nrow = 2) %>% round(digits = 3)
##        [,1]  [,2]
## [1,]  0.936 0.289
## [2,] -0.223 0.182
opt_mod_sd_correctlyspecified$B %>% matrix(nrow = 2) %>% round(digits = 3)
##       [,1]  [,2]
## [1,] 0.693 1.088
## [2,] 0.279 0.260
opt_mod_sd_misspecified$B %>% matrix(nrow = 2) %>% round(digits = 3)
##       [,1]  [,2]
## [1,] 0.638 0.376
## [2,] 0.287 0.176

Print distribution parameters with non-robust and robust standard errors , \(\lambda, p, q\), for each variable (in columns)

opt_mod_params$distr %>% matrix(ncol = 2) %>% round(digits = 3)
##           [,1]  [,2]
## [1,]    -0.524 0.149
## [2,]     1.492 1.921
## [3,] 85616.311 7.891
opt_mod_sd_correctlyspecified$distr %>% matrix(ncol = 2) %>% round(digits = 3)
##       [,1]    [,2]
## [1,] 1.094   1.419
## [2,] 2.826  12.576
## [3,] 0.000 409.497
opt_mod_sd_misspecified$distr %>% matrix(ncol = 2) %>% round(digits = 3)
##       [,1]    [,2]
## [1,] 0.444   2.295
## [2,] 2.261  21.203
## [3,] 0.000 614.852

5.1 Alternative Option

Under correct specification may also choose, e.g., the analytic version of the Hessian to obtain the information matrix and calculate a robust version by using again the analytic version of the Hessian as well as the OPG

opt_mod_params = fill_tmpl_whf_rev(theta = PARAMS,
                                   tmpl = TMPL)
opt_mod_sd_correctlyspecified_alt = fill_tmpl_whf_rev(theta = comp_sd %>% pull(ac_ha_inv),
                                                  tmpl = TMPL)
opt_mod_sd_misspecified_alt = fill_tmpl_whf_rev(theta = comp_sd %>% pull(bb_rob_ha),
                                            tmpl = TMPL)

Print AR parameters with non-robust and robust standard errors in wide matrix format.

opt_mod_params$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]   [,3]   [,4]
## [1,]    1    0 -0.107 -0.277
## [2,]    0    1  0.492 -1.013
opt_mod_sd_correctlyspecified_alt$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 2.626 1.294
## [2,]    0    1 0.847 0.289
opt_mod_sd_misspecified_alt$polm_ar %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 4.801 2.473
## [2,]    0    1 1.095 0.187

Print MA backwards parameters with non-robust and robust standard errors in wide matrix format.

opt_mod_params$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 0.065 0.749
## [2,]    0    1 0.259 0.002
opt_mod_sd_correctlyspecified_alt$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 2.332 1.691
## [2,]    0    1 0.967 1.604
opt_mod_sd_misspecified_alt$polm_ma_bwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 3.933 1.605
## [2,]    0    1 1.409 3.492

Print MA forwards parameters with non-robust and robust standard errors in wide matrix format.

opt_mod_params$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]   [,3]   [,4]
## [1,]    1    0  0.561  1.671
## [2,]    0    1 -0.241 -0.296
opt_mod_sd_correctlyspecified_alt$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 1.055 3.111
## [2,]    0    1 0.353 1.272
opt_mod_sd_misspecified_alt$polm_ma_fwd %>% matrix(nrow = 2) %>% round(digits = 3)
##      [,1] [,2]  [,3]  [,4]
## [1,]    1    0 1.588 4.675
## [2,]    0    1 0.510 2.671

Print static matrix B with non-robust and robust standard errors.

opt_mod_params$B %>% matrix(nrow = 2) %>% round(digits = 3)
##        [,1]  [,2]
## [1,]  0.936 0.289
## [2,] -0.223 0.182
opt_mod_sd_correctlyspecified_alt$B %>% matrix(nrow = 2) %>% round(digits = 3)
##       [,1]  [,2]
## [1,] 0.644 0.616
## [2,] 0.258 0.198
opt_mod_sd_misspecified_alt$B %>% matrix(nrow = 2) %>% round(digits = 3)
##       [,1]  [,2]
## [1,] 0.644 0.364
## [2,] 0.279 0.177

Print distribution parameters with non-robust and robust standard errors , \(\lambda, p, q\), for each variable (in columns)

opt_mod_params$distr %>% matrix(ncol = 2) %>% round(digits = 3)
##           [,1]  [,2]
## [1,]    -0.524 0.149
## [2,]     1.492 1.921
## [3,] 85616.311 7.891
opt_mod_sd_correctlyspecified_alt$distr %>% matrix(ncol = 2) %>% round(digits = 3)
##       [,1]    [,2]
## [1,] 0.683   1.690
## [2,] 2.379  15.377
## [3,] 0.000 467.316
opt_mod_sd_misspecified_alt$distr %>% matrix(ncol = 2) %>% round(digits = 3)
##       [,1]    [,2]
## [1,] 0.469   2.295
## [2,] 2.124  21.204
## [3,] 0.000 614.855

6 Skewness and Kurtosis

sgt_skewness = function(l, p, q){
  v = q^(-1/p) *
    (
      (3*l^2+1) *
        (beta(3/p,q-2/p)/beta(1/p,q)) -
        4*l^2 *
        (beta(2/p,q-1/p)/beta(1/p,q))^2
    )^(-1/2)
  
  skewness = 2 * q^(3/p) * l * v^3 / beta(1/p, q)^3 *
    (8*l^2*beta(2/p,q-1/p)^3 - 
       3*(1+3*l^2) * beta(1/p,q) * beta(2/p, q-1/p) * beta(3/p,q-2/p) + 
       2 * (1+l^2) * beta(1/p,q)^2 * beta(4/p, q-3/p)
       )
  
  return(skewness)
}

sgt_kurtosis = function(l, p, q){
  v = q^(-1/p) *
    (
      (3*l^2+1) *
        (beta(3/p,q-2/p)/beta(1/p,q)) -
        4*l^2 *
        (beta(2/p,q-1/p)/beta(1/p,q))^2
    )^(-1/2)
  
  kurtosis = q^(4/p) * v^4 / beta(1/p, q)^4 *
    (48*l^4*beta(2/p,q-1/p)^4 +
       24 * l^2 * (1+3*l^2) * beta(1/p,q) * beta(2/p, q-1/p)^2 * beta(3/p,q-2/p) -
       32 * l^2 * (1+l^2) * beta(1/p,q)^2 * beta(2/p, q-1/p) * beta(4/p,q-3/p) + 
       (1+10*l^2+5*l^4) * beta(1/p,q)^3 * beta(5/p, q-4/p)
    )
  
  return(kurtosis)
}

sgt_skewness(-0.52,1.49,85616)
## [1] -1.035358
sgt_kurtosis(-0.52,1.49,85616)
## [1] 6.318257

sgt_skewness(0.15,1.92,7.89)
## [1] 0.3251307
sgt_kurtosis(0.15,1.92,7.89)
## [1] 3.769704